///////////////////////////////////////////////////////////////////////////////////////////////////
/// \file formic/fqmc/vmc_walker.cpp
///
/// \brief   implementation file for the variational monte carlo walker class
///
///////////////////////////////////////////////////////////////////////////////////////////////////

#include<string>
#include<vector>
#include<complex>
#include<cmath>
#include<algorithm>
#include<numeric>

#include<formic/fqmc/vmc_walker.h>
#include<formic/wfn/symmetry.h>
#include<formic/input/base.h>
#include<formic/numeric/numeric.h>
#include<formic/output/output.h>
#include<formic/timing/timing.h>
#include<formic/fqmc/distance_proposer.h>
#include<formic/mpi/interface.h>
#include<formic/wfn/configuration.h>

///////////////////////////////////////////////////////////////////////////////////////////////////
/// \brief Constructs the walker.
///
/// \param[in]      userinp  the user's input options
///
///////////////////////////////////////////////////////////////////////////////////////////////////
template <class S> formic::VMC_Walker<S>::VMC_Walker(const formic::InputBase & userinp)
  : _accepted(false),
    _count(1.0),
    _nacc(0.0),
    _npro(0.0),
    _mps(userinp.get<int>("moves_per_sample")),
    _steps(0),
    _nirr(userinp.get<std::vector<int> >("opirr").size()),
    _config(userinp),
    _config_diff(),
    _dpt(),
    _orb_irrep(userinp.get<std::vector<int> >("orb_irrep")),
    _lfg(),
    _gf(),
    _co_vec(),
    _accumulators(),
    _dist_prop(userinp)
{

  // get the point group's direct product table
  formic::group_direct_product_table(userinp.get<std::string>("point_group"), _dpt);

//  // get the irrep for each orbital
//  {
//    const std::vector<int> & opirr = userinp.get<std::vector<int> >("opirr");
//    int p = 0;
//    for (int i = 0; i < opirr.size(); i++)
//      for (int j = 0; j < opirr.at(i); j++)
//        _orb_irrep.at(p++) = i;
//    if ( p != _orb_irrep.size() )
//      throw formic::Exception("counting error when assigning orbital irreps");
//  }

}

template formic::VMC_Walker<double>::VMC_Walker(const formic::InputBase &);
template formic::VMC_Walker<std::complex<double> >::VMC_Walker(const formic::InputBase &);

///////////////////////////////////////////////////////////////////////////////////////////////////
/// \brief  Proposes a configuration difference corresponding to a possible walker move
///
/// \param[in]      userinp  the user's input options
///
/// \return  the proposal density ratio
///
///////////////////////////////////////////////////////////////////////////////////////////////////
template <class S> double formic::VMC_Walker<S>::propose_config_diff(const formic::InputBase & userinp) {

  // get number of atomic orbitals
  const int nao = userinp.get<int>("nao");

  // get the number of alpha and beta electrons
  const int na = userinp.get<int>("nalpha");
  const int nb = userinp.get<int>("nbeta");

  // get references to the occ/unocc maps
  const std::vector<int> & ao2s = _config.a_occ_to_site();
  const std::vector<int> & bo2s = _config.b_occ_to_site();
  const std::vector<int> & au2s = _config.a_unocc_to_site();
  const std::vector<int> & bu2s = _config.b_unocc_to_site();

  bool bad_move = false;

  // construct moves at random until a suitable move is found
  do {

    bad_move = false;

    // pick two different electrons at random
    int elec1 = _lfg() % ( na + nb - 0 );
    int elec2 = _lfg() % ( na + nb - 1 );
    elec2 += ( elec2 >= elec1 ? 1 : 0 );

    // get the spins of the chosen electrons
    int spin1 = ( elec1 < na ? 0 : 1 );
    int spin2 = ( elec2 < na ? 0 : 1 );

    // convert the electron indices now that we know their spin
    elec1 -= spin1 * na;
    elec2 -= spin2 * na;

    // get the lattice sites the electrons occupy
    int pos1 = ( spin1 == 0 ? 2*ao2s.at(elec1)+0 : 2*bo2s.at(elec1)+1 );
    int pos2 = ( spin2 == 0 ? 2*ao2s.at(elec2)+0 : 2*bo2s.at(elec2)+1 );
    if ( pos2 < pos1) {
      std::swap( pos1,  pos2);
      std::swap(elec1, elec2);
      std::swap(spin1, spin2);
    }

    // get an integer used when the spins are the same
    int same_spin = ( spin1 == spin2 ? 1 : 0 );

    // pick two holes at random to move the electrons into
    int hole1 = _lfg() % ( ( spin1 == 0 ? au2s.size() : bu2s.size() ) - 0         );
    int hole2 = _lfg() % ( ( spin2 == 0 ? au2s.size() : bu2s.size() ) - same_spin );
    hole2 += ( same_spin && hole2 >= hole1 ? 1 : 0 );

    // get the new positions
    int new1 = ( spin1 == 0 ? 2*au2s.at(hole1)+0 : 2*bu2s.at(hole1)+1 );
    int new2 = ( spin2 == 0 ? 2*au2s.at(hole2)+0 : 2*bu2s.at(hole2)+1 );

    // check that point group symmetry is satisfied
    int irrep1 = _dpt.at( _orb_irrep.at(pos1/2) + _nirr * _orb_irrep.at(new1/2) );
    int irrep2 = _dpt.at( _orb_irrep.at(pos2/2) + _nirr * _orb_irrep.at(new2/2) );
    bad_move = bad_move || irrep1 != irrep2;

    // define the difference object
    _config_diff.newly_empty().resize(2);
    _config_diff.newly_empty().at(0) = pos1;
    _config_diff.newly_empty().at(1) = pos2;
    _config_diff.newly_full().resize(2);
    _config_diff.newly_full().at(0) = new1;
    _config_diff.newly_full().at(1) = new2;

    //if (!bad_move)
    //  formic::of << boost::format("moving %2i  %2i    to    %2i  %2i") % pos1 % pos2 % new1 % new2
    //             << std::endl
    //             << std::endl;

  } while (bad_move);

  // for now, the proposal density is symmetric, so return one
  return 1.0;

}

///////////////////////////////////////////////////////////////////////////////////////////////////
/// \brief Proposes a move to a new configuration, which is assumed to be stored internally.
///
/// \param[in]      userinp  the user's input options
///
///////////////////////////////////////////////////////////////////////////////////////////////////
template <class S> void formic::VMC_Walker<S>::propose_move(const formic::InputBase & userinp) {

  // propose where to move to and get the proposal density ratio
//  const double pdr = formic::Configuration::propose_move(userinp, _lfg, _config, _config_diff);
//  const double pdr = this->propose_config_diff(userinp);
  const double pdr = _dist_prop.propose(userinp, _config, _config_diff, _lfg);

  // calculate the ratio of the guiding function's weight for the new and old configurations
  S gf_ratio = _gf->get_ratio(_config_diff);

  //const S gf_new_val = gf_ratio * _gf->get_value().get();
  //const S gf_old_val = _gf->get_value().get();

  //const formic::RestPairingFunc<S> * const rpf = dynamic_cast<formic::RestPairingFunc<S> *>( _gf.get() );
  //if ( !rpf->wfn_sampling() )
  //  gf_ratio = formic::unity<S>() * std::sqrt( std::sqrt( formic::square_norm(gf_ratio) ) );

  // choose a random number and decide whether to accept the proposed move according to the Metropolis algorithm
  _accepted = ( pdr != 0.0 && formic::uni_01(_lfg) < pdr * formic::square_norm(gf_ratio) );

  // update propose and accept counts
  _npro += 1.0;
  if ( _accepted )
    _nacc += 1.0;

//        // check the value of the old configuration
//        if ( false && !rpf->wfn_sampling() ) {
//          const S old_val = _gf->get_value().get();
//          const S cor_val = rpf->correct_value(_config, 'v');
//          //std::cout << boost::format("    proposed  %20.12f          correct  %20.12f           %s")
//          //             % old_val
//          //             % cor_val
//          //             % ( _accepted ? "accepted" : "" )
//          //          << std::endl;
//          assert( std::abs( old_val - cor_val ) < 1.0e-9 );
//        }
//
//        // check the value of the new configuration
//        if ( false && !rpf->wfn_sampling() && _accepted ) {
//          const S cor_val = rpf->correct_value(_next_config, 'v');
//          std::cout << boost::format("  proposed  %20.12f     correct  %20.12f           %s")
//                       % gf_new_val
//                       % cor_val
//                       % ( _accepted ? "accepted" : "" )
//                    ;
//                    //<< std::endl;
//          assert( std::abs( gf_new_val - cor_val ) < 1.0e-9 );
//        }
//
//        if (true) {
//          //int unpaired = 0;
//          //for (int i = 0; i < _config.fock_config().size(); i+=2)
//          //  if ( _config.fock_config().at(i) != _config.fock_config().at(i+1) )
//          //    unpaired++;
//          //formic::of << "unpaired/2 odd ?     " << ((unpaired/2) % 2) << std::endl;
//          for (int i = 0; i < _config.fock_config().size(); i+=2)
//            formic::of << _config.fock_config().at(i);
//          formic::of << " ";
//          for (int i = 1; i < _config.fock_config().size(); i+=2)
//            formic::of << _config.fock_config().at(i);
//          //formic::of << boost::format("      newly occ:  %2i  %2i   newly emp:  %2i  %2i  ")
//          //              % _config_diff.newly_full().at(0)
//          //              % _config_diff.newly_full().at(1)
//          //              % _config_diff.newly_empty().at(0)
//          //              % _config_diff.newly_empty().at(1);
//          formic::of << boost::format("    pdr  %20.12f          square wfn ratio  %20.12f           %s")
//                       % pdr
//                       % formic::square_norm(gf_ratio)
//                       % ( _accepted ? "accepted" : "" )
//                    << std::endl;
//        }

}

template void formic::VMC_Walker<double>::propose_move(const formic::InputBase &);
template void formic::VMC_Walker<std::complex<double> >::propose_move(const formic::InputBase &);

///////////////////////////////////////////////////////////////////////////////////////////////////
/// \brief Accumulate data at the current configuration.
///
/// \param[in]      userinp  the user's input options
///
///////////////////////////////////////////////////////////////////////////////////////////////////
template <class S> void formic::VMC_Walker<S>::accumulate(const formic::InputBase & userinp) {

  //if ( _steps != 0 )
  //  return;

  // if the move was accepted or we are doing more than one move per sample, accumulate data for the configuration the walker is about to leave
  if ( _mps > 1 || _accepted ) {
    typename std::vector<boost::shared_ptr<formic::Accumulator<S> > >::iterator it = _accumulators.begin();
    for ( ; it != _accumulators.end(); it++)
      (**it).accumulate(userinp, _count);
  }

}

template void formic::VMC_Walker<double>::accumulate(const formic::InputBase &);
template void formic::VMC_Walker<std::complex<double> >::accumulate(const formic::InputBase &);

///////////////////////////////////////////////////////////////////////////////////////////////////
/// \brief Move the walker according to the proposed move and whether it was accepted.
///
/// \param[in]      userinp  the user's input options
///
///////////////////////////////////////////////////////////////////////////////////////////////////
template <class S> void formic::VMC_Walker<S>::move(const formic::InputBase & userinp) {

  //formic::start_timer("VMC_Walker::move");

  // if the move was accepted...
  if (_accepted) {

    // update the guiding function for the new configuration
    _gf->update(_config_diff);

    //std::cout << boost::format("move accepted,  new coeff = %20.12f        noiseless coeff = %20.12f")
    //             % _gf->get_value().get()
    //             % _gf->noiseless_value(_next_config).get()
    //          << std::endl
    //          << std::endl;

    // update the non-guiding coefficient value objects for the new configuration
    typename std::vector<boost::shared_ptr<formic::CoeffObj<S> > >::iterator it = _co_vec.begin();
    for ( ; it != _co_vec.end(); it++)
      if (*it != _gf)
        (**it).update(_config_diff);

    // update the current configuration to reflect the move
    _config.update(_config_diff);

    // reset the counter for how long we have been on the current configuration
    _count = 1.0;

  // if the move was rejected...
  } else {

    //std::cout << "move rejected" << std::endl << std::endl;

    // update the counter for how long we have been on the current configuration
    _count += 1.0;

  }

  //_steps = ( _steps + 1 ) % userinp.nparticles();
  //_count = 1.0;

  // the count should always be one when making more than one move per sample
  if ( _mps > 1 )
    _count = 1.0;

  //formic::stop_timer("VMC_Walker::move");

}

template void formic::VMC_Walker<double>::move(const formic::InputBase &);
template void formic::VMC_Walker<std::complex<double> >::move(const formic::InputBase &);

///////////////////////////////////////////////////////////////////////////////////////////////////
/// \brief Perform any remaining accumulations and communicate the accumulators' results.
///
/// \param[in]      userinp  the user's input options
///
///////////////////////////////////////////////////////////////////////////////////////////////////
template<class S> void formic::VMC_Walker<S>::final_accumulate(const formic::InputBase & userinp) {

  // accumulate data for the final configuration
  const double count_offset = 1.0; //( _mps > 1 ? 0.0 : 1.0 );
  for (typename std::vector<boost::shared_ptr<formic::Accumulator<S> > >::iterator it = _accumulators.begin();
       it != _accumulators.end();
       it++)
    (**it).accumulate(userinp, _count - count_offset);
    //(**it).accumulate(userinp, _count);

  // finilize the accumulators
  for (typename std::vector<boost::shared_ptr<formic::Accumulator<S> > >::iterator it = _accumulators.begin();
       it != _accumulators.end();
       it++)
    (**it).finalize(userinp);

  // get and print overall acceptance rate
  const int nred = 2;
  double x[nred];
  double y[nred];
  x[0] = _nacc;
  x[1] = _npro;
  formic::mpi::reduce(x, y, nred, MPI::SUM);
  if ( formic::mpi::rank() == 0 )
    formic::of << boost::format("sampling finished.    acceptance ratio = %.12f") % ( y[0] / y[1] ) << std::endl;

}

template void formic::VMC_Walker<double>::final_accumulate(const formic::InputBase &);
template void formic::VMC_Walker<std::complex<double> >::final_accumulate(const formic::InputBase &);

///////////////////////////////////////////////////////////////////////////////////////////////////
/// \brief reset the walker in preparation for sampling
///
/// \param[in]      userinp  the user's input options
///
///////////////////////////////////////////////////////////////////////////////////////////////////
template <class S> void formic::VMC_Walker<S>::reset(const formic::InputBase & userinp, const bool reset_co) {

  // reset coefficient objects if requested
  if (reset_co) {

    // initialize the guiding function
    _gf->initialize(_config);

    // initialize the non-guiding coefficient value objects
    for (typename std::vector<boost::shared_ptr<formic::CoeffObj<S> > >::iterator it = _co_vec.begin();
         it != _co_vec.end();
         it++)
      if (*it != _gf)
        (**it).initialize(_config);

  }

  // reset accept and propose counts
  _nacc = 0.0;
  _npro = 0.0;

  // reset the accumulators
  for (typename std::vector<boost::shared_ptr<formic::Accumulator<S> > >::iterator it = _accumulators.begin();
       it != _accumulators.end();
       it++)
    (**it).reset(_config);
  //struct local_func {
  //  const formic::Configuration & _cfg;
  //  local_func(const formic::Configuration & config) : _cfg(config) {}
  //  void operator()(boost::shared_ptr<formic::Accumulator<S> > & acc) { acc->reset(_cfg); }
  //};
  //std::for_each(_accumulators.begin(), _accumulators.end(), local_func(_config));

  // initialize the count of how long we have been at the current configuration
  _count = 1.0;

  //// reset the number of steps taken since the last accumulation
  //_steps = 0;

}

template void formic::VMC_Walker<double>::reset(const formic::InputBase &, const bool);
template void formic::VMC_Walker<std::complex<double> >::reset(const formic::InputBase &, const bool);

///////////////////////////////////////////////////////////////////////////////////////////////////
/// \brief  take an exact sample using the guiding function and accumulators
///
/// \param[in]      userinp  the user's input options
///
///////////////////////////////////////////////////////////////////////////////////////////////////
template<class S> void formic::VMC_Walker<S>::exact_sample(const formic::InputBase & userinp) {

  // if necessary, prepare exact sample configuration index lists
  if ( this->gf()->exact_ia_list().size() == 0 || this->gf()->exact_ib_list().size() == 0 )
    this->gf()->prep_for_exact_sample(userinp);

  // create a flag to tell when to reset the accumulators
  bool have_reset_acc = false;

  //int to_send = 0;
  //if ( formic::mpi::size() == 2 && formic::mpi::rank() != 0 )
  //  formic::mpi::recv(&to_send, 1, 0, 0);

  const bool detailed_config_print = userinp.get<bool>("print_config_info");
  if ( detailed_config_print && formic::mpi::size() > 1 )
    throw formic::Exception("printing config info in exact sample only works for mpi::size() == 1");

  // loop over this process's configurations
  std::vector<int>::const_iterator ia = this->gf()->exact_ia_list().begin();
  std::vector<int>::const_iterator fa = this->gf()->exact_ia_list().end();
  std::vector<int>::const_iterator ib = this->gf()->exact_ib_list().begin();
  for ( ; ia != fa; ia++, ib++) {

    // populate alpha configuration
    this->config().populate_alpha_or_beta(0, *ia);

    // populate beta configuration
    this->config().populate_alpha_or_beta(1, *ib);

    formic::start_timer("exact sample prep coeff obj");

    // initialize coefficient objects for this configuration
    this->gf()->initialize(this->config());
    for (int i = 0; i < this->co_vec().size(); i++)
      if ( this->co_vec().at(i) != this->gf() )
        this->co_vec().at(i)->initialize(this->config());

    // get the coefficient value for the configuration
    const S cv = this->gf()->get_coeff().get();

    formic::stop_timer("exact sample prep coeff obj");

    // if we have not already reset the accumulators, do so now
    if ( !have_reset_acc ) {
      for (int i = 0; i < this->accumulators().size(); i++)
        this->accumulators().at(i)->reset(this->config());
      have_reset_acc = true;
    }

    // set the weight equal to the square of the coefficient value
    _count = formic::square_norm(cv);

    if ( detailed_config_print )
      formic::of << boost::format("  %14.4e  %14.4e  ") % formic::square_norm(cv) % cv;

    // accumulate for this configuration
    _accepted = true;
    this->accumulate(userinp);

    // print configuration
    if ( detailed_config_print ) {
      const int no = this->config().fock_config().size()/2;
      //formic::of << "    ";
      //for(int i = 0; i < no; i++)
      //  formic::of << "  " << this->config().fock_config().at(2*i+0) << this->config().fock_config().at(2*i+1);
      //formic::of << boost::format("  %20.12f") % cv << std::endl;
      //formic::of.flush();
      formic::of << "    ";
      for(int i = 0; i < no; i++)
        formic::of << "  " << this->config()[2*i+0] << this->config()[2*i+1];
      formic::of << "\n";
      //formic::of.flush();
    }

  }

  if ( detailed_config_print )
    formic::of << std::endl;

  //if ( formic::mpi::size() == 2 && formic::mpi::rank() == 0 )
  //  formic::mpi::send(&to_send, 1, 1, 0);

  // finalize the accumulators
  for (int i = 0; i < this->accumulators().size(); i++)
    this->accumulators().at(i)->finalize(userinp);

//  // local functions
//  struct local_funcs {
//
//    // function to populate the alpha or beta configuration based on the compound index x
//    static void populate_config(const int no, const int b, int x, formic::Configuration & config) {
//      for (int i = 0; i < no; i++) {
//        config[2*i+b] = x - 2 * ( x / 2 ); //x % 2;
//        x /= 2;
//      }
//    }
//
//    // counts the number of alpha or beta electrons in a configuration
//    static int count_pop(const int no, const int b, const formic::Configuration & config) {
//      int retval = 0;
//      for (int i = 0; i < no; i++) {
//        assert( config[2*i+b] == 0 || config[2*i+b] == 1 );
//        retval += config[2*i+b];
//      }
//      return retval;
//    }
//
//    // returns the irreducible representation of the configuration
//    static int config_irrep(const formic::Configuration & config, const int nirr, const std::vector<int> & orb_irrep, const std::vector<int> & dir_prod_table) {
//      int retval = 0;
//      for (int i = 0; i < config.fock_config().size(); i++)
//        retval = ( config[i] == 0 ? retval : dir_prod_table.at( orb_irrep.at(i/2) + nirr * retval ) );
//      return retval;
//    }
//
//  };
//
//  // get the number of orbitals
//  const int no = userinp.get<int>("nao");
//
//  // get the number of alpha/beta electrons
//  const int na = userinp.get<int>("nalpha");
//  const int nb = userinp.get<int>("nbeta");
//
//  // get the number of irreps
//  const int nirr = userinp.get<std::vector<int> >("opirr").size();
//
//  // get the orbital irreps
//  const std::vector<int> & orb_irrep = userinp.get<std::vector<int> >("orb_irrep");
//
//  // get the point group's direct product table
//  std::vector<int> dir_prod_table;
//  formic::fqmc::group_direct_product_table(userinp.get<std::string>("point_group"), dir_prod_table);
//
//  // get the irrep of the initial configuration, which sets the wavefunction's irrep
//  const int target_irrep = local_funcs::config_irrep(formic::Configuration(userinp), nirr, orb_irrep, dir_prod_table);
//
//  // get the maximum spin configuration compound index
//  const int maxi = formic::pow<int>(2, no);
//
//  // get the small coefficient cutoff
//  const double coeff_cutoff = userinp.get<double>("exact_sample_cutoff");
//
//  // First, count the non-zero configurations on the root node
//  int nnonzero = 0;
//  if (formic::mpi::rank() == 0) {
//
//    // initialize a counter for the number of configurations skipped due to small coefficients
//    int small_skip_count = 0;
//
//    // vector of coefficient norms above the cutoff
//    std::vector<double> coeff_norms;
//
//    // loop over alpha configurations
//    for (int ia = 0; ia < maxi; ia++) {
//
//      // populate the configuration
//      local_funcs::populate_config(no, 0, ia, this->config());
//
//      //formic::of << "    " << boost::format("%10i") % ia;
//      //formic::of << "    ";
//      //for(int i = 0; i < no; i++)
//      //  formic::of << this->config()[2*i+0];
//      //formic::of.flush();
//
//      // skip configurations with the wrong number of alpha particles
//      if ( local_funcs::count_pop(no, 0, this->config()) != na ) {
//        //formic::of << "  skipped" << std::endl;
//        continue;
//      }
//
//      //formic::of << std::endl;
//
//      // loop over beta configurations
//      for (int ib = 0; ib < maxi; ib++) {
//
//        // populate the configuration
//        local_funcs::populate_config(no, 1, ib, this->config());
//
//        //// print configuration
//        //formic::of << "    " << boost::format("%10i") % ia;
//        //formic::of << "    ";
//        //for(int i = 0; i < no; i++)
//        //  formic::of << this->config()[2*i+0];
//        //formic::of << "    " << boost::format("%10i") % ib;
//        //formic::of << "    ";
//        //for(int i = 0; i < no; i++)
//        //  formic::of << this->config()[2*i+1];
//        //formic::of.flush();
//
//        // skip configurations with the wrong number of beta particles
//        if ( local_funcs::count_pop(no, 1, this->config()) != nb ) {
//          //formic::of << "  skipped" << std::endl;
//          continue;
//        }
//
//        // skip configurations with the wrong irrep
//        if ( target_irrep != local_funcs::config_irrep(this->config(), nirr, orb_irrep, dir_prod_table) ) {
//          //formic::of << "  skipped" << std::endl;
//          continue;
//        }
//
//        //// print configuration
//        ////formic::of << "    " << boost::format("%10i") % ia;
//        //formic::of << "    ";
//        //for(int i = 0; i < no; i++)
//        //  formic::of << this->config()[2*i+0];
//        ////formic::of << "    " << boost::format("%10i") % ib;
//        //formic::of << "    ";
//        //for(int i = 0; i < no; i++)
//        //  formic::of << this->config()[2*i+1];
//        //formic::of.flush();
//
//        formic::start_timer("exact sample prep coeff obj");
//
//        // initialize the configuration
//        this->config().init_occ_unocc_site();
//
//        // get the coefficient value for the configuration
//        const S cv = this->gf()->get_coeff(this->config()).get();
//
//        //// print coefficient value
//        //formic::of << "    " << formic::format_number("%20.12f", cv) << std::endl;
//
//        formic::stop_timer("exact sample prep coeff obj");
//
//        // skip configurations with zero coefficients
//        if ( std::abs(cv) < coeff_cutoff ) {
//          small_skip_count++;
//          //for(int i = 0; i < no; i++)
//          //  formic::of << this->config()[2*i+0];
//          ////formic::of << "    " << boost::format("%10i") % ib;
//          //formic::of << "    ";
//          //for(int i = 0; i < no; i++)
//          //  formic::of << this->config()[2*i+1];
//          //formic::of << "    " << formic::format_number("%20.6e", cv) << "     skipped" << std::endl;
//          continue;
//        }
//
//        // count the configuration if it's coefficient is nonzero
//        nnonzero++;
//
//        // record the coefficient norm
//        coeff_norms.push_back( std::abs(cv) );
//
//      }
//
//    }
//
//    // sort the coefficient norms by size (making the median easy to find)
//    std::sort(coeff_norms.begin(), coeff_norms.end());
//
//    // print some info on the sample
//    formic::of << boost::format("skipped %i configurations due to small coefficients during exact sampling") % small_skip_count << std::endl;
//    formic::of << boost::format("  info on norms above cutoff:") << std::endl;
//    formic::of << boost::format("         minimum = %.2e") % *coeff_norms.begin() << std::endl;
//    formic::of << boost::format("         maximum = %.2e") % *coeff_norms.rbegin() << std::endl;
//    formic::of << boost::format("         average = %.2e") % ( std::accumulate(coeff_norms.begin(), coeff_norms.end(), 0.0) / coeff_norms.size() ) << std::endl;
//    formic::of << boost::format("          median = %.2e") % coeff_norms.at(coeff_norms.size()/2) << std::endl;
//
//  }
//
//  // send the number of nonzero configurations to all nodes
//  formic::mpi::bcast(nnonzero);
//
//  // create a flag to tell when to reset the accumulators
//  bool have_reset_acc = false;
//
//  // determine which configurations to do on this process
//  const int np = formic::mpi::size();
//  const int mr = formic::mpi::rank();
//  const int my_first_config = ( nnonzero / np ) * mr;
//  const int my_last_config = ( mr < np - 1 ? ( nnonzero / np ) * ( mr + 1 ) : nnonzero ) - 1;
//
//  // create a counter for the number of nonzero configurations seen
//  int count_nonzero = 0;
//
//  // loop over alpha configurations
//  for (int ia = 0; ia < maxi; ia++) {
//
//    // populate the configuration
//    local_funcs::populate_config(no, 0, ia, this->config());
//
//    //formic::of << "    " << boost::format("%10i") % ia;
//    //formic::of << "    ";
//    //for(int i = 0; i < no; i++)
//    //  formic::of << this->config()[2*i+0];
//    //formic::of.flush();
//
//    // skip configurations with the wrong number of alpha particles
//    if ( local_funcs::count_pop(no, 0, this->config()) != na ) {
//      //formic::of << "  skipped" << std::endl;
//      continue;
//    }
//
//    //formic::of << std::endl;
//
//    // loop over beta configurations
//    for (int ib = 0; ib < maxi; ib++) {
//
//      // populate the configuration
//      local_funcs::populate_config(no, 1, ib, this->config());
//
//      //// print configuration
//      //formic::of << "    " << boost::format("%10i") % ia;
//      //formic::of << "    ";
//      //for(int i = 0; i < no; i++)
//      //  formic::of << this->config()[2*i+0];
//      //formic::of << "    " << boost::format("%10i") % ib;
//      //formic::of << "    ";
//      //for(int i = 0; i < no; i++)
//      //  formic::of << this->config()[2*i+1];
//      //formic::of.flush();
//
//      // skip configurations with the wrong number of beta particles
//      if ( local_funcs::count_pop(no, 1, this->config()) != nb ) {
//        //formic::of << "  skipped" << std::endl;
//        continue;
//      }
//
//      // skip configurations with the wrong irrep
//      if ( target_irrep != local_funcs::config_irrep(this->config(), nirr, orb_irrep, dir_prod_table) ) {
//        //formic::of << "  skipped" << std::endl;
//        continue;
//      }
//
//      //// print configuration
//      ////formic::of << "    " << boost::format("%10i") % ia;
//      //formic::of << "    ";
//      //for(int i = 0; i < no; i++)
//      //  formic::of << this->config()[2*i+0];
//      ////formic::of << "    " << boost::format("%10i") % ib;
//      //formic::of << "    ";
//      //for(int i = 0; i < no; i++)
//      //  formic::of << this->config()[2*i+1];
//      //formic::of.flush();
//
//      formic::start_timer("exact sample prep coeff obj");
//
//      // initialize the configuration
//      this->config().init_occ_unocc_site();
//
//      // get the coefficient value for the configuration
//      const S cv = this->gf()->get_coeff(this->config()).get();
//
//      //// print coefficient value
//      //formic::of << "    " << formic::format_number("%20.12f", cv) << std::endl;
//
//      // skip configurations with zero coefficients
//      if ( std::abs(cv) < coeff_cutoff ) {
//        formic::stop_timer("exact sample prep coeff obj");
//        //for(int i = 0; i < no; i++)
//        //  formic::of << this->config()[2*i+0];
//        ////formic::of << "    " << boost::format("%10i") % ib;
//        //formic::of << "    ";
//        //for(int i = 0; i < no; i++)
//        //  formic::of << this->config()[2*i+1];
//        //formic::of << "    " << formic::format_number("%20.6e", cv) << "     skipped" << std::endl;
//        continue;
//      }
//
//      // only do this process's configurations
//      if ( count_nonzero < my_first_config || count_nonzero > my_last_config ) {
//        formic::stop_timer("exact sample prep coeff obj");
//        count_nonzero++;
//        continue;
//      }
//      count_nonzero++;
//
//      // initialize coefficient objects for this configuration
//      this->gf()->initialize(this->config());
//      for (int i = 0; i < this->co_vec().size(); i++)
//        if ( this->co_vec().at(i) != this->gf() )
//          this->co_vec().at(i)->initialize(this->config());
//
//      formic::stop_timer("exact sample prep coeff obj");
//
//      // if we have not already reset the accumulators, do so now
//      if ( !have_reset_acc ) {
//        for (int i = 0; i < this->accumulators().size(); i++)
//          this->accumulators().at(i)->reset(this->config());
//        have_reset_acc = true;
//      }
//
//      //// print configuration
//      ////formic::of << "after reset:" << std::endl;
//      ////formic::of << "    " << boost::format("%10i") % ia;
//      ////formic::of << "    ";
//      //for(int i = 0; i < no; i++)
//      //  formic::of << this->config()[2*i+0];
//      ////formic::of << "    " << boost::format("%10i") % ib;
//      //formic::of << "    ";
//      //for(int i = 0; i < no; i++)
//      //  formic::of << this->config()[2*i+1];
//      //formic::of << "    " << formic::format_number("%20.6e", cv) << std::endl;
//      ////formic::of << std::endl;
//
//      // set the weight equal to the square of the coefficient value
//      _count = formic::square_norm(cv);
//
//      // accumulate for this configuration
//      _accepted = true;
//      this->accumulate(userinp);
//
//      //// print configuration
//      //formic::of << "after accumulate:" << std::endl;
//      //formic::of << "    " << boost::format("%10i") % ia;
//      //formic::of << "    ";
//      //for(int i = 0; i < no; i++)
//      //  formic::of << this->config()[2*i+0];
//      //formic::of << "    " << boost::format("%10i") % ib;
//      //formic::of << "    ";
//      //for(int i = 0; i < no; i++)
//      //  formic::of << this->config()[2*i+1];
//      //formic::of << std::endl;
//
//    }
//
//  }
//
//  //formic::of << std::endl;
//
//  // finalize the accumulators
//  for (int i = 0; i < this->accumulators().size(); i++)
//    this->accumulators().at(i)->finalize(userinp);

}

template void formic::VMC_Walker<double>::exact_sample(const formic::InputBase &);
template void formic::VMC_Walker<std::complex<double> >::exact_sample(const formic::InputBase &);
